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We present some preliminary results from recent numerical simulations that model the evolution of super-massive black 
hole (SMBH) binaries in galactic nuclei. Including the post-Newtonian terms for the binary system and adopting appro- 
priate models for the galaxies allows us, for the first time, to follow the evolution of SMBH binaries from kpc scales down 
to the coalescence phase. We use our results to make predictions of the detectability of such events with the gravitational 
wave detector LISA. 
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1 Introduction 

The formation and evolution of super-massive black holes 
(SMBHs) in galactic nuclei is an important ingredient for 
our understanding of galaxy formation and evolution. In the 
hierarchical cosmological framework, galaxies merge in the 
course of their evolution. Therefore, either many massive 
binary (or even multiple) SMBHs should exist, or else we 
need to develop a detailed understanding of what ultimately 
happens to the massive binaries. According to the standard 
theory, the evolution of SMBHs after a galaxy merger can 
be divided in three stages (Begelman, Blandford & Rees 
1980): (i) Dynamical friction causes the individual SMBHs 
to spiral to the galactic centre where they form a binary, (ii) 
Super-elastic scattering of stars causes the binary to harden, 
(iii) Eventually, the binary may coalesce owing to the emis- 
sion of gravitational waves. However, this is only possible if 
stellar- and gas-dynamical processes have first brought the 
two SMBHs to small enough separations (some 10 _3 pc) 
that the emission of gravitational waves is significant. Wheth-| 
er Nature typically succeeds in overcoming this "final par- 
sec problem," (see review by Milosavljevic & Merritt 2003) 
is currently unknown. 

Here we present preliminary results of recent A-body 
simulations of SMBH binaries in galactic nuclei, including 
post-Newtonian CPA) equations of motion up to the 2.57W 
order. In Sec. 2 we describe our VAf implementation and 
provide numerical tests and simple two-body experiments. 
In Sec. 3 we summarise some first results of our detailed 
galactic nuclei simulations, following the evolution of the 



SMBH binaries from the unbound state down to relativistic 
coalescence. We conclude with Sec. 4. 



2 A two-body Hermite VM integrator 

In order to incorporate the dynamics of relativistic binary 
systems into our direct A-body code (see Sec. 3) we use the 
"PA^equations of motion, which are expressed by expanding 
the relativistic acceleration between two compact objects in 
a power series of 1/c (e.g., Damour & Deruelle 1981; Sof- 
fel 1989). Schematically, the VM equation of motion for an 
object in the binary can be written as: 
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where ao is the classical Newtonian acceleration, ai and 
&2 are the energy conserving TPA/'and 27W'corrections, re- 
spectively, and a 2 .5 is the dissipative 2.57-Wterm, describ- 
ing the dominant damping forces due to gravitational wave 
emission. In the current VM routine we neglect any effects 
of spin-orbit and spin-spin coupling (e.g., Damour 1987). 
For the full expressions of the individual terms up to 3.5VJV 
order see, e.g., Blanchet (2006). We implement the equa- 
tions of motion formulated in a general harmonic (Carte- 
sian) coordinate system, whereas other astrophysical imple- 
mentations often use expressions given in the binary's cen- 
tre of mass frame (e.g., Kupi, Amaro-Seoane & Spurzem 
2006; Aarseth 2007). 

In order to achieve high accuracy in the simulations we 
use a 4 th order Hermite integrator (e.g., Makino & Aarseth 
1992), which requires in addition to the accelerations also 
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Fig. 1 Comparison of the VM implementation (solid line) and 
the finite differences Acn/Ati (symbols). Shown are the curves 
for the 27Wx-component of some (arbitrary) orbit integration. 
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Fig. 2 Inspiral and plunge of an equal-mass binary system using 
the post-Newtonian equations of motion. The numerical integra- 
tion is done using a 4 th order Hermite integrator and hierarchical 
time-step scheme. 

their first time-derivatives ("jerk") - including the VM cor- 
rections. Since the VM accelerations and jerks are compli- 
cated functions of the particles' masses, positions and ve- 
locities, their derivation and numerical implementation is a 
potential source of error. Therefore, we provide some nu- 
merical tests which turned out to be useful benchmarks, as 
already applied with other TWimplementations (e.g., Kupi 
et al. 2006; Lockmann & Baumgardt 2008): 

(i) Time-derivatives. To check the correct implementa- 
tion of the jerks of the "PA^terms, we evolve a binary system 
for several orbital periods and compare the finite differences 
Aa-pj^/At for consecutive time-steps to the directly calcu- 
lated a from our implementation (Fig. 1). Any deviations 
between the two would be easily detectable and clearly in- 
dicate any errors in the values of the computed jerk. In Fig. 2 
we show a typical trajectory of an binary inspiral integrated 
with our TWcode, starting with a quasi-circular orbit. 

(ii) Merging time for eccentric binary orbits. In the 
following test we set the gravitational constant G and the 
two particle masses to unity. The speed of light in our model 
units is set to c = 10. 1 The two particles are placed on the 
x-axis with an initial separation of Ax = 2. The initial ve- 
locities ±v y are calculated for a Keplerian orbit of given ec- 



These units are chosen for numerical testing purposes only and are not 
physically motivated. 



Table 1 The columns from left to right give for each integrated 
orbit the corresponding (1 — e%), the initial eccentricity e , the 
initial semi-major axis ao, and the two 'merging' times TvM and 
T2.5VM, respectively. 
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centricity e . We then integrate the orbits using (a) the full 
VM (i.e., in our case up to 2.5VJV) corrections and (b) ap- 
plying only the 2.57'A^correction. In both cases the integra- 
tions are stopped at the times T-p^f and T2.5VM, respectively, 
when the binaries reach a separation of 10 Schwarzschild 
radii (Rg = 2Gmbh/c 2 ), i.e., well before the VM approx- 
imation may become inadequate. The results for different 
initial eccentricities are given in Table 1 . The numbers given 
in Table 1 are in agreement with integrations using the VM 
implementations of Kupi et al. (2006) and Aarseth (2007). 
The exact result may slightly vary for different integration 
schemes and time-step controls. 

In Fig. 3 we show a comparison of the "merging" times 
T-ptv (solid line) and T 2 . (dashed line). Note that the re- 
sults can vary by a factor of a few with increasing eccentric- 
ity. This underlines the importance of the lower order VN 
corrections, which sometimes have been neglected in ear- 
lier astrophysical simulations. The pericentre shift due to 
the WNanA 27^ terms leads to more pericentre passages 
and therefore to an overall stronger dissipation (especially 
for high eccentricities), which shortens T-pj^f as compared 
to T 2 .hVN'- The almost horizontal part of the two curves in 
Fig. 3 (1 — eg « 0.2. . .0.4) is the result of a direct plunge into 
the 10 i?s region due to the high eccentricity and small peri- 
centre distances. Finally, we compare our results to the nu- 
merical integration of the equations given by Peters (1964), 
describing the orbit-averaged changes of the eccentricity 
and semi-major axis due to gravitational wave emission. We 
find good agreement between the latter integration and our 
corresponding 2.5VN two-body integrations. The diverg- 
ing behaviour for (1 — eg) < 0.4 comes from the fact that 
our orbit-averaged integration is terminated when the semi- 
major axis (instead of the binary separation) reaches 10 i?s- 

(iii) The innermost stable circular orbit (ISCO). Here 
we study the linear stability of circular orbits. It is known 
for binary systems that circular orbits with radii below a 
certain radius n are unstable, even without dissipative ef- 
fects. The evolution of an relativistic binary changes from a 
quasi-adiabatic inspiral to an unstable plunge when crossing 
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Fig. 3 Shown are the results obtained from numerical integration 
of the orbit-averaged equations from Peters (1964) (dotted line) 
and from the direct two-body integration including the 2.57Wcor- 
rection (dashed line) and full 'PA/'corrections. 

n (e.g., Kidder, Will & Wiseman 1993). Analytically, the 
radius n of the ISCO can be calculated using the standard 
linear stability analysis (e.g., Blanchet & Iyer 2003). To de- 
termine n on an ISCO for objects of comparable mass we 
use here the exact 2VAf accurate circular frequencies u^vM 
as given in Preto (2007). We find for the binary system de- 
scribed in test (ii) that the ISCO has a radius of r\ ~ 0. 136 in 
our model units (Preto 2007). Now we determine r\ by nu- 
merical integration of circular VM orbits. The model units 
and masses are the same as the ones used in test (ii). The 
initial velocities of the two particles are calculated based 
on the TPM accurate u)rpH (see above). We integrate the 
VM circular orbits including only YPM and 2VM correc- 
tions, i.e., under the absence of gravitational wave emission. 
The integration of the orbits has been terminated at a time 
*torm either after an unstable plunge or after a total maxi- 
mum of rii (initial) orbital periods T orb . In Fig. 4 we show 
ttcrm/irii x T or b) for rii — 100 (solid line), rii = 500 (dot- 
ted line) and rii = 1000 (dashed line). As an indicator for 
stability we also plot the final ratio r tcim /r , where r tcrm 
is the separation of the particles at time ttorm- We find that 
our VM two-body integration accurately captures the radius 
of the ISCO in all three cases. The transition from the stable 
to the unstable regime is very abrupt as the theory predicts, 
and happens at the radius estimated by the linear stability 
analysis. 

3 Super-massive black hole binaries in 
galactic nuclei 

In this section we present preliminary results of recent direct 
A^-body simulations of SMBH binaries in galactic nuclei in- 
cluding VM corrections (Berentzen et al. 2008). Initial con- 
ditions were based on the models of Berczik et al. (2006). 
The initial field particle distribution was a discrete reali- 
sation of the axisymmetric, rotating King (1966) models 
of Longaretti & Lagoute (1996). We consider such mod- 



Fig. 4 Normalized integration time t tclnl /(ni x T or b) for ini- 
tially circular orbits with radius ro. We show our results for n» = 
100 (solid line), m = 500 (dotted line) and n t = 1000 (dashed 
line). The symbols show the final ratio r term /ro. The vertical line 
indicates the radius of the innermost stable circular orbit (ISCO) 
as derived from the linear stability analysis. Note that our VMcode 
is capable of accurately capturing the ISCO. 

els to be reasonable representations of galactic nuclei in 
post-merger galaxies. When the rotation parameter is suf- 
ficiently large, the models evolve rapidly into triaxial bars. 
Berczik et al. (2006) found that in this case, the massive bi- 
nary continued to harden at an essentially A^-independent 
rate, due presumably to the centrophilic nature of the or- 
bits (Merritt & Poon 2004). This feature means that evo- 
lution of the massive binary can be faithfully reproduced 
even with modest particle numbers. Accordingly, the nu- 
cleus was represented here using either A^ = 25 x 10 3 or 
N = 50 x 10 3 particles, in each case with nine different 
random realisations. The masses of the two SMBH parti- 
cles were chosen to be one percent of the total mass. For 
the simulations we used an updated version of the publi- 
cally available 2 ip -GRAPE code (Harfst et al. 2007). The 
code is an NBODYl-like algorithm (Aarseth 1999), includ- 
ing a hierarchical time-step scheme and a 4* h -order Her- 
mite integrator. Our modified version includes the VM im- 
plementation described in the previous section. Our simu- 
lations were carried out on the high-performance GRAPE- 
6A clusters at the Astronomisches Rechen-Institut (Heidel- 
berg), Rochester Institute of Technology (New York) and 
the Main Astron. Observatory (Kiev). More details on the 
code and hardware are given in Berentzen et al. (2008) and 
Spurzem et al. (2008). 

In Fig. 5 we show the evolution of the inverse semi- 
major axis 1/a and eccentricity e of the binaries in our 
simulations with full 7W(solid lines) and with only 2.5VM 
corrections (dashed lines). Initially the two SMBH particles 
are unbound; they form a very eccentric binary at t ~ 25. 
Afterwards the binaries steadily harden mainly by stellar- 
dynamical processes, before relativistic effects start to dom- 
inate their evolution. Eventually, they reach the inspiral pha- 

2 http://wiki.cs.rit.edu/view/GRAPEcluster/phiGRAPE 
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Fig. 5 Evolution of the SMBH binaries in a model of a rotating 
galactic nucleus. Shown are the orbital eccentricity e (thin lines) 
and the inverse semi-major axis (thick lines) as a function of time. 
The full and dashed lines represent models using the full VM and 
2. 5VAf corrections, respectively. Note the significantly different 
evolution of the binary. 
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Fig. 6 Merging time r morg as a function of eccentricity e. Shown 
are simulations with 25k (circles) and 50k (squares) particles. The 
open and filled symbols represent models with full VM and only 
2. 57W corrections, respectively. The gray region marks the theo- 
retically estimated merging times (see Berentzen et al. 2008). 

se during which the orbits circularise before the final coales- 
cence. As expected, we find that the lower order VM terms 
significantly influence the evolution of the binaries and must 
not be neglected in this type of simulation. Finally, we show 
in Fig. 6 the merging times T mer g for our SMBH binaries as 
a function of e. Scaling our model units to real galaxies we 
find that r merg in our models is typically less than 1 Gyr 
after the binary has formed. 

4 Summary 

We upgraded the direct iV-body code tp-GRAPE by includ- 
ing a relativistic treatment for binary SMBHs via the VN 
equations of motion up to 2.57Worder. The new implemen- 
tation was carefully tested for accuracy. We then carried out 
numerical simulations of SMBH binaries in galactic nuclei. 
Our models allow us, for the first time, to follow the evolu- 



tion of initially unbound SMBHs with kiloparsec-scale sep- 
arations down to relativistic inspiral and coalescence. The 
SMBH binaries in our simulations form generally with high 
eccentricities, which are maintained during the evolution- 
ary phase dominated by stellar-dynamical effects. We have 
shown that the inclusion of the lower order YPNwd TPM 
corrections is crucial for obtaining the correct time depen- 
dence of the binary orbital parameters and final coalescence 
time. Our results demonstrate that SMBH binaries in certain 
models of galactic nuclei can overcome the stalling barrier 
by stellar-dynamical effects alone. A detailed description 
and discussion of our models and results will be given in 
Berentzen et al. (2008) and in Preto et al. (2008). 
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